
module sulchem

   !----------------------------------------------------------------------- 
   ! Purpose: 
   ! Contains the sulfur cycle chemistry codes.
   ! 
   ! Author: 
   ! Chemistry code is from Mary Barth.
   ! Module coded by B. Eaton.
   !----------------------------------------------------------------------- 

   use shr_kind_mod, only: r8 => shr_kind_r8
   use ppgrid,       only: pcols, pver
   use abortutils,   only: endrun
   use physconst,    only: tmelt
   use perf_mod
   use cam_logfile,  only: iulog

   implicit none
   save
   private
   public :: &
      inisulchem,   &! Initialize sulfur cycle chemistry.
      sulf_chemdr,  &
      cldychmini, dicor

   ! private module data

   integer, parameter ::&
      nz=18,              &
      nza=8,              &
      nalb=2,             &
      nz1=nz-1,           &
      nza1=nza-1,         &
      nz2=51,             &
      nsiz =nz *nza*nalb, &
      nsiz2=nz2*nza*nalb

   real(r8) ::&
      jper2(nsiz2),          &! j values interpolated to new ht grid
      delz(nz1),             &! Difference between tabl. z's
      delang(nza1),          &! Difference between tabl. sec(zen.ang)
      delalb,                &! Difference in tabl. albedo = .5-.2
      pie,                   &! 3.14159...
      dayspy                  ! Number of days per 1 year

   real(r8), parameter :: mwa2so2 = 28.966_r8/64._r8 ! ratio molecular wt of air to so2
   real(r8), parameter :: mwa2h2o2 = 28.966_r8/34._r8 ! ratio molecular wt of air to h2o2


   real(r8), parameter, dimension(nsiz) :: jper =    &! Tabulated j values
(/-11.63010025_r8,-11.53059959_r8,-11.47570038_r8,-11.44130039_r8,-11.40859985_r8,-11.40690041_r8,-11.42399979_r8,-11.46020031_r8,&
  -11.48509979_r8,-11.49100018_r8,-11.46759987_r8,-11.41129971_r8,-11.32050037_r8,-11.20450020_r8,-10.96660042_r8,-10.65079975_r8,&
  -10.31389999_r8,-10.15680027_r8,-11.93789959_r8,-11.81239986_r8,-11.73900032_r8,-11.68900013_r8,-11.62989998_r8,-11.60620022_r8,&
  -11.60490036_r8,-11.61820030_r8,-11.62650013_r8,-11.62040043_r8,-11.58790016_r8,-11.52509975_r8,-11.43029976_r8,-11.31110001_r8,&
  -11.06599998_r8,-10.75179958_r8,-10.39190006_r8,-10.20240021_r8,-12.20409966_r8,-12.06019974_r8,-11.97229958_r8,-11.90939999_r8,&
  -11.82750034_r8,-11.78409958_r8,-11.76560020_r8,-11.75730038_r8,-11.74950027_r8,-11.73139954_r8,-11.68959999_r8,-11.62010002_r8,&
  -11.52070045_r8,-11.39859962_r8,-11.14780045_r8,-10.83170033_r8,-10.46100044_r8,-10.24470043_r8,-12.51459980_r8,-12.35309982_r8,&
  -12.24969959_r8,-12.17240047_r8,-12.06379986_r8,-11.99619961_r8,-11.95629978_r8,-11.92070007_r8,-11.89220047_r8,-11.85890007_r8,&
  -11.80539989_r8,-11.72710037_r8,-11.62139988_r8,-11.49520016_r8,-11.23849964_r8,-10.91720009_r8,-10.54150009_r8,-10.29699993_r8,&
  -13.14570045_r8,-12.95919991_r8,-12.83010006_r8,-12.72649956_r8,-12.56420040_r8,-12.44439983_r8,-12.35589981_r8,-12.25739956_r8,&
  -12.18089962_r8,-12.11240005_r8,-12.03289986_r8,-11.93480015_r8,-11.81389999_r8,-11.67720032_r8,-11.40919971_r8,-11.07330036_r8,&
  -10.69830036_r8,-10.41100025_r8,-14.34010029_r8,-14.15550041_r8,-14.01509953_r8,-13.88949966_r8,-13.65480042_r8,-13.43579960_r8,&
  -13.23670006_r8,-12.97910023_r8,-12.77519989_r8,-12.61489964_r8,-12.47119999_r8,-12.32730007_r8,-12.17119980_r8,-12.00699997_r8,&
  -11.70969963_r8,-11.35060024_r8,-10.96660042_r8,-10.65600014_r8,-15.41930008_r8,-15.26990032_r8,-15.16590023_r8,-15.07689953_r8,&
  -14.89920044_r8,-14.68239975_r8,-14.41380024_r8,-13.96300030_r8,-13.55609989_r8,-13.23820019_r8,-12.98649979_r8,-12.77019978_r8,&
  -12.56400013_r8,-12.36279964_r8,-12.01889992_r8,-11.63490009_r8,-11.22630024_r8,-10.92059994_r8,-17.23080063_r8,-17.09090042_r8,&
  -16.99939919_r8,-16.92959976_r8,-16.82430077_r8,-16.74169922_r8,-16.66230011_r8,-16.49410057_r8,-16.09469986_r8,-15.37609959_r8,&
  -14.65410042_r8,-14.07880020_r8,-13.62609959_r8,-13.26119995_r8,-12.75020027_r8,-12.25749969_r8,-11.79109955_r8,-11.45209980_r8,&
  -11.17109966_r8,-11.14850044_r8,-11.14080048_r8,-11.14089966_r8,-11.15730000_r8,-11.18900013_r8,-11.23019981_r8,-11.29240036_r8,&
  -11.33609962_r8,-11.35649967_r8,-11.34700012_r8,-11.30399990_r8,-11.22620010_r8,-11.12189960_r8,-10.90159988_r8,-10.60239983_r8,&
  -10.27840042_r8,-10.12619972_r8,-11.52420044_r8,-11.47599983_r8,-11.44919968_r8,-11.43290043_r8,-11.42099953_r8,-11.42879963_r8,&
  -11.44970036_r8,-11.48649979_r8,-11.51080036_r8,-11.51659966_r8,-11.49520016_r8,-11.44289970_r8,-11.35809994_r8,-11.24800014_r8,&
  -11.01669979_r8,-10.71510029_r8,-10.36550045_r8,-10.18050003_r8,-11.81890011_r8,-11.75230026_r8,-11.71059990_r8,-11.68099976_r8,&
  -11.64519978_r8,-11.63220024_r8,-11.63479996_r8,-11.64830017_r8,-11.65499973_r8,-11.64729977_r8,-11.61489964_r8,-11.55410004_r8,&
  -11.46290016_r8,-11.34809971_r8,-11.10849953_r8,-10.80249977_r8,-10.44029999_r8,-10.22789955_r8,-12.15250015_r8,-12.06820011_r8,&
  -12.01109982_r8,-11.96700001_r8,-11.90400028_r8,-11.86610031_r8,-11.84650040_r8,-11.83150005_r8,-11.81620026_r8,-11.79199982_r8,&
  -11.74639988_r8,-11.67520046_r8,-11.57610035_r8,-11.45580006_r8,-11.20800018_r8,-10.89459991_r8,-10.52560043_r8,-10.28439999_r8,&
  -12.80440044_r8,-12.69659996_r8,-12.61540031_r8,-12.54609966_r8,-12.43089962_r8,-12.34109974_r8,-12.27270031_r8,-12.19359970_r8,&
  -12.12870026_r8,-12.06789970_r8,-11.99429989_r8,-11.90139961_r8,-11.78509998_r8,-11.65240002_r8,-11.39009953_r8,-11.05939960_r8,&
  -10.68850040_r8,-10.40359974_r8,-13.98359966_r8,-13.87979984_r8,-13.79139996_r8,-13.70489979_r8,-13.52700043_r8,-13.34549999_r8,&
  -13.17129993_r8,-12.93599987_r8,-12.74400043_r8,-12.59039974_r8,-12.45110035_r8,-12.31060028_r8,-12.15719986_r8,-11.99520016_r8,&
  -11.70090008_r8,-11.34430027_r8,-10.96230030_r8,-10.65270042_r8,-15.04539967_r8,-14.97089958_r8,-14.91469955_r8,-14.86209965_r8,&
  -14.74139977_r8,-14.57019997_r8,-14.33740044_r8,-13.92109966_r8,-13.53120041_r8,-13.22140026_r8,-12.97410011_r8,-12.76060009_r8,&
  -12.55630016_r8,-12.35649967_r8,-12.01439953_r8,-11.63179970_r8,-11.22420025_r8,-10.91909981_r8,-16.85479927_r8,-16.78750038_r8,&
  -16.74150085_r8,-16.70509911_r8,-16.64769936_r8,-16.59880066_r8,-16.54520035_r8,-16.40889931_r8,-16.04310036_r8,-15.35249996_r8,&
  -14.64319992_r8,-14.07279968_r8,-13.62240028_r8,-13.25860023_r8,-12.74860001_r8,-12.25650024_r8,-11.79049969_r8,-11.45160007_r8 /)

   real(r8), parameter, dimension(nz) :: zz =     &! Tabulated heights
      (/ 0._r8, 1._r8, 2._r8, 3._r8, 5._r8, 7._r8, 9._r8, 12._r8, 15._r8, 18._r8, 21._r8, 24._r8, &
         27._r8, 30._r8, 35._r8, 40._r8, 45._r8, 50._r8 /)

   real(r8), parameter, dimension(nza) :: zasec = &! Tabulated sec(zen.ang)
      (/ 1._r8, 1.3_r8, 1.6_r8, 2._r8, 3._r8, 6._r8, 12._r8, 50._r8 /)

   real(r8), parameter, dimension(nalb) :: alb =  &! Tabulated albedo = 0.2 and 0.5
      (/ 0.2_r8, 0.5_r8 /)

   real(r8), parameter, dimension(121) :: yield = &! yield of DMS+OH rxn going to SO2 
      (/ 0.76905_r8,  0.76953_r8,  0.77002_r8,  0.77050_r8,  0.77098_r8,  0.77146_r8, &
         0.77195_r8,  0.77243_r8,  0.77291_r8,  0.77340_r8,  0.77388_r8,  0.77436_r8, &
         0.77485_r8,  0.77533_r8,  0.77581_r8,  0.77629_r8,  0.77678_r8,  0.77726_r8, &
         0.77774_r8,  0.77822_r8,  0.77870_r8,  0.77918_r8,  0.77967_r8,  0.78015_r8, &
         0.78063_r8,  0.78111_r8,  0.78159_r8,  0.78207_r8,  0.78255_r8,  0.78303_r8, &
         0.78351_r8,  0.78399_r8,  0.78447_r8,  0.78495_r8,  0.78543_r8,  0.78591_r8, &
         0.78640_r8,  0.78688_r8,  0.78737_r8,  0.78786_r8,  0.78835_r8,  0.78884_r8, &
         0.78934_r8,  0.78984_r8,  0.79035_r8,  0.79086_r8,  0.79137_r8,  0.79190_r8, &
         0.79243_r8,  0.79297_r8,  0.79353_r8,  0.79409_r8,  0.79467_r8,  0.79526_r8, &
         0.79587_r8,  0.79650_r8,  0.79715_r8,  0.79783_r8,  0.79853_r8,  0.79927_r8, &
         0.80004_r8,  0.80084_r8,  0.80169_r8,  0.80258_r8,  0.80352_r8,  0.80452_r8, &
         0.80558_r8,  0.80671_r8,  0.80790_r8,  0.80918_r8,  0.81055_r8,  0.81200_r8, &
         0.81356_r8,  0.81522_r8,  0.81700_r8,  0.81890_r8,  0.82093_r8,  0.82309_r8, &
         0.82540_r8,  0.82786_r8,  0.83047_r8,  0.83325_r8,  0.83618_r8,  0.83928_r8, &
         0.84255_r8,  0.84598_r8,  0.84957_r8,  0.85332_r8,  0.85722_r8,  0.86126_r8, &
         0.86543_r8,  0.86973_r8,  0.87412_r8,  0.87861_r8,  0.88316_r8,  0.88777_r8, &
         0.89240_r8,  0.89705_r8,  0.90169_r8,  0.90631_r8,  0.91087_r8,  0.91537_r8, &
         0.91978_r8,  0.92409_r8,  0.92829_r8,  0.93236_r8,  0.93630_r8,  0.94009_r8, &
         0.94373_r8,  0.94721_r8,  0.95054_r8,  0.95370_r8,  0.95670_r8,  0.95954_r8, &
         0.96223_r8,  0.96476_r8,  0.96714_r8,  0.96938_r8,  0.97148_r8,  0.97344_r8, &
         0.97528_r8 /)

   real(r8), parameter, dimension(121) :: dmsrate = &! rate coeff of DMS+OH 
      (/ 0.21261E-10_r8,0.21112E-10_r8,0.20966E-10_r8,0.20823E-10_r8,0.20683E-10_r8, &
         0.20546E-10_r8,0.20411E-10_r8,0.20280E-10_r8,0.20151E-10_r8,0.20024E-10_r8, &
         0.19900E-10_r8,0.19779E-10_r8,0.19660E-10_r8,0.19543E-10_r8,0.19428E-10_r8, &
         0.19316E-10_r8,0.19206E-10_r8,0.19098E-10_r8,0.18991E-10_r8,0.18887E-10_r8, &
         0.18785E-10_r8,0.18684E-10_r8,0.18585E-10_r8,0.18488E-10_r8,0.18393E-10_r8, &
         0.18299E-10_r8,0.18207E-10_r8,0.18116E-10_r8,0.18027E-10_r8,0.17939E-10_r8, &
         0.17852E-10_r8,0.17766E-10_r8,0.17682E-10_r8,0.17598E-10_r8,0.17516E-10_r8, &
         0.17434E-10_r8,0.17354E-10_r8,0.17274E-10_r8,0.17194E-10_r8,0.17115E-10_r8, &
         0.17036E-10_r8,0.16958E-10_r8,0.16880E-10_r8,0.16801E-10_r8,0.16722E-10_r8, &
         0.16643E-10_r8,0.16563E-10_r8,0.16482E-10_r8,0.16401E-10_r8,0.16317E-10_r8, &
         0.16233E-10_r8,0.16146E-10_r8,0.16057E-10_r8,0.15966E-10_r8,0.15871E-10_r8, &
         0.15774E-10_r8,0.15673E-10_r8,0.15567E-10_r8,0.15458E-10_r8,0.15343E-10_r8, &
         0.15223E-10_r8,0.15097E-10_r8,0.14965E-10_r8,0.14826E-10_r8,0.14680E-10_r8, &
         0.14526E-10_r8,0.14365E-10_r8,0.14195E-10_r8,0.14016E-10_r8,0.13828E-10_r8, &
         0.13632E-10_r8,0.13426E-10_r8,0.13211E-10_r8,0.12987E-10_r8,0.12754E-10_r8, &
         0.12514E-10_r8,0.12265E-10_r8,0.12009E-10_r8,0.11747E-10_r8,0.11479E-10_r8, &
         0.11207E-10_r8,0.10932E-10_r8,0.10655E-10_r8,0.10376E-10_r8,0.10098E-10_r8, &
         0.98217E-11_r8,0.95482E-11_r8,0.92787E-11_r8,0.90145E-11_r8,0.87565E-11_r8, &
         0.85057E-11_r8,0.82630E-11_r8,0.80289E-11_r8,0.78042E-11_r8,0.75893E-11_r8, &
         0.73845E-11_r8,0.71899E-11_r8,0.70058E-11_r8,0.68321E-11_r8,0.66687E-11_r8, &
         0.65155E-11_r8,0.63722E-11_r8,0.62384E-11_r8,0.61140E-11_r8,0.59985E-11_r8, &
         0.58915E-11_r8,0.57926E-11_r8,0.57013E-11_r8,0.56173E-11_r8,0.55402E-11_r8, &
         0.54694E-11_r8,0.54047E-11_r8,0.53456E-11_r8,0.52917E-11_r8,0.52426E-11_r8, &
         0.51981E-11_r8,0.51578E-11_r8,0.51214E-11_r8,0.50886E-11_r8,0.50591E-11_r8, &
         0.50327E-11_r8 /)


!##############################################################################
contains
!##############################################################################


  subroutine inisulchem()

    use cam_history,    only: addfld, add_default, phys_decomp

    !----------------------------------------------------------------------- 
    ! Purpose: 
    ! Initialize sulfur cycle chemistry module.
    ! 
    ! Author: B. Eaton
    !----------------------------------------------------------------------- 

    implicit none

    ! Local variables
    integer ::&
         i, l, ik, k

    !-----------------------------------------------------------------------
    !     initialize variables for jh2o2 interpolation
    !     interpolate jval data onto 1 km height grid
    do l = 1, nza*nalb
       do ik = 1, nz2-1
          k = 2
11        continue
          if( zz(k) .gt. ik-1 ) then
             jper2((l-1)*nz2+ik) = jper((l-1)*nz+k-1) &
                  + ((ik-1)-zz(k-1))/(zz(k)-zz(k-1)) &
                  * (jper((l-1)*nz+k)-jper((l-1)*nz+k-1))
          else
             k = k + 1
             go to 11
          endif
       end do
       jper2(l*nz2) = jper(l*nz) 
    end do

    do i = 1,nza1
       delang(i) = 1._r8 / (zasec(i+1) - zasec(i))
    end do
    delalb  = 1._r8 / (alb(2) - alb(1))

    !     Constants from radiation routines needed for dicor()
    dayspy  =  365._r8
    pie     =  4._r8*atan(1._r8)

    call addfld('DMSSNK  ','kg/kg/s' ,pver, 'A', 'DMSsink ',phys_decomp)
    call add_default ('DMSSNK', 1, ' ')
    call addfld('SO2SRCG ','kg/kg/s' ,pver, 'A', 'SO2 Src G   ',phys_decomp)
    call add_default ('SO2SRCG', 1, ' ')
    call addfld('SO2SRCG2','kg/kg/s' ,pver, 'A', 'SO2 Src G2  ',phys_decomp)
    call add_default ('SO2SRCG2', 1, ' ')
    call addfld('SO4SRC  ','kg/kg/s' ,pver, 'A', 'SO4 Src Tot ',phys_decomp)
    call add_default ('SO4SRC', 1, ' ')
    call addfld('SO4SRCG ','kg/kg/s' ,pver, 'A', 'SO4 Src G   ',phys_decomp)
    call add_default ('SO4SRCG', 1, ' ')
    call addfld('H2O2SRC ','kg/kg/s' ,pver, 'A', 'H2O2Src ',phys_decomp)
    call add_default ('H2O2SRC', 1, ' ')
    call addfld('H2O2SNKA','kg/kg/s' ,pver, 'A', 'H2O2SnkA',phys_decomp)
    call add_default ('H2O2SNKA', 1, ' ')
    call addfld('H2O2SNKG','kg/kg/s' ,pver, 'A', 'H2O2SnkG',phys_decomp)
    call add_default ('H2O2SNKG', 1, ' ')

    call addfld ('PH','pH',pver, 'A','Cloud Water pH',phys_decomp)
    call add_default ('PH', 1, ' ')

  end subroutine inisulchem

!##############################################################################
   subroutine sulf_chemdr( t, pmid, cldwat, rain, cldf, cldv, &
                           icwmr1, icwmr2, indcp, ncldypts, hion, &
                           so2, so4, dms,  h2o2, ekso2, ekh2o2, &
                           o3, h2o23d, oh, no3, ho2, q,  &
                           zm, dtime, calday, clat, clon, ncol, lchnk  )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Sulfur chemistry driver separated from chemwdepdr.  The prognostic species 
      ! are DMS, SO2, SO4, and H2O2.  Needed a chem driver without the wet 
      ! deposition for MOZART chemistry
      ! 
      ! Author: F Vitt
      !----------------------------------------------------------------------- 

      use cam_history, only: outfld

      real(r8), intent(in) ::&
         clat(pcols),              &
         clon(pcols),              &
         calday,                   &
         dtime,                    &! time step
         pmid(pcols,pver),         &! pressure at layer midpoints
         zm(pcols,pver),           &! height of layer midpoints
         t(pcols,pver),            &! air temperature (K)
         q(pcols,pver),            &! specific humidity (kg/kg)
         cldwat(pcols,pver),       &! cloud water mixing ratio
         cldf(pcols,pver),         &! total cloud fraction
         cldv(pcols,pver),         &! cloudy volume which is occupied by rain or cloud water 
         rain(pcols,pver),         &! total precip mixing ratio
         icwmr1(pcols,pver),       &! in cloud water mixing ratio for zhang+hack scheme
         icwmr2(pcols,pver)         ! in cloud water mixing ratio for zhang+hack scheme

      real(r8), dimension(pcols,pver), intent(inout) ::&
         so2,         &! so2 array
         so4,         &! so4 array
         dms,         &! dms
         h2o2          ! h2o2 

      real(r8), dimension(pcols,pver), intent(in) ::&
         o3,          &
         no3,         &
         ho2,         &
         oh,          &
         h2o23d

      integer, intent(in) ::    ncol
      integer, intent(in) ::    lchnk

      real(r8), intent(in) ::&
         ekh2o2(pcols,pver)    ! H2O2  Henry's Law coefficient
      real(r8), intent(inout) ::&
         ekso2(pcols,pver)     ! SO2  effective Henry's Law coefficient

      integer, intent(in) ::&
         indcp(pcols,pver),    &! Indices of cloudy grid points
         ncldypts(pver)         ! Number of cloudy grid points

      real(r8), dimension(pcols,pver), intent(inout) ::&
         hion        ! Hydrogen ion concentration in drops

      ! Local variables:

      real(r8), dimension(pcols,pver) ::&
         h2o2new,     &! h2o2
         h2o2tmp,     &! temporary h2o2
         so2new,      &! so2 array
         so2tmp,      &! temporary so2 array
         so4new,      &! so4 array
         so4tmp,      &! temporary so4 array
         dmsnew,      &! dms
         dmssnk,      &! total sink of dms
         so2srcg,     &! total source of so2  
         so2srcg2,    &! dms + oh source of so2  
         so4src,      &! total source of so4  
         so4srcg,     &! gas source of so4  
         so4srca,     &! aq. source of so4  
         h2o2src,     &! total source of h2o2  
         h2o2snkg,    &! gas sink of h2o2  
         h2o2snka      ! aqueous sink of h2o2  

      real(r8), dimension(pcols,pver) ::&
         jh2o2         ! photolysis rate of H2O2

      integer :: i, k

!----------------------------------------------------------------------

     call getj( clat, clon, calday, zm, ncol, jh2o2 )

     call aqchem( so2, so4, h2o2, t, pmid, &
          cldwat, rain, so2new, so4new, h2o2new, &
          o3, cldf, so4srca, hion, h2o2snka, &
          ekso2, ekh2o2, dtime, indcp, ncldypts, &
          cldv, icwmr1, icwmr2, ncol )

     so2tmp(:ncol,:) = so2new(:ncol,:)
     so4tmp(:ncol,:) = so4new(:ncol,:)
     h2o2tmp(:ncol,:) = h2o2new(:ncol,:)

     call gaschem( so2tmp, so4tmp, dms, h2o2tmp, t, &
          pmid, so2new, so4new, dmsnew, h2o2new, &
          h2o23d, oh, no3, ho2, q, &
          jh2o2, so4srcg, dmssnk, so2srcg, so2srcg2, &
          h2o2src, h2o2snkg, dtime, ncol )

#ifdef SCYC_MASSBGT
         call endrun ('sulf_chemdr: gamwet expects scalar for actual argument lat')
         call gamwet( 'dmssnk', lat, pdel, dmssnk )
         call gamwet( 'so2srcg', lat, pdel, so2srcg )
         call gamwet( 'so4srca', lat, pdel, so4srca )
         call gamwet( 'so4srcg', lat, pdel, so4srcg )
#endif

      so2( :ncol,:) = so2new( :ncol,:)
      so4( :ncol,:) = so4new( :ncol,:)
      dms( :ncol,:) = dmsnew( :ncol,:)
      h2o2(:ncol,:) = h2o2new(:ncol,:)

      do k = 1, pver
         do i = 1, ncol
            so4src(i,k) = so4srca(i,k) + so4srcg(i,k)
         end do
      end do

      call outfld( 'DMSSNK  ', dmssnk, pcols, lchnk )
      call outfld( 'SO2SRCG ', so2srcg, pcols, lchnk )
      call outfld( 'SO2SRCG2', so2srcg2, pcols, lchnk )
      call outfld( 'SO4SRC  ', so4src, pcols, lchnk )
      call outfld( 'SO4SRCG ', so4srcg, pcols, lchnk )
      call outfld( 'H2O2SRC ', h2o2src, pcols, lchnk )
      call outfld( 'H2O2SNKA', h2o2snka, pcols, lchnk )
      call outfld( 'H2O2SNKG', h2o2snkg, pcols, lchnk )

   endsubroutine sulf_chemdr


   subroutine aqchem( so2, so4, h2o2, tair, pres, &
                      qc, qr, so2new, so4new, h2o2new, &
                      o3, cldf, so4srca, hion, h2o2snka, &
                      ekso2, ekh2o2, ztodt, ind, ncpts, &
                      cldv, icwmr1, icwmr2, ncol )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Calculate the sources and sinks from the aqueous chemistry:
      !          S(IV) + H2O2 --> SO4 
      !          S(IV) + O3   --> SO4 
      !
      !  where concentrations are concentrations in aqueous phase.
      !
      ! Author: M. Barth
      !----------------------------------------------------------------------- 

      implicit none

      integer, intent(in) ::&
!         lat(pcols),          &! latitude indices
!         lon(pcols),          &! longititude indices
         ind(pcols,pver),     &! indices for cloudy grid points
         ncpts(pver),         &! number of cloudy grid points
         ncol                  ! number of atmospheric columns used in chunk
      real(r8), intent(in) ::&
         so2(pcols,pver),     &! so2 array
         so4(pcols,pver),     &! so4 from aqueous chemistry 
         h2o2(pcols,pver),    &! h2o2 in kg/kg
         o3(pcols,pver),      &! o3
         tair(pcols,pver),    &! air temperature (K)
         pres(pcols,pver),    &! midpoint pressures
         qc(pcols,pver),      &! cloud water mixing ratio
         qr(pcols,pver),      &! rain mixing ratio
         ztodt,               &! time step
         cldf(pcols,pver)     ! fraction cloud cover

      real(r8), intent(in) ::&
         cldv(pcols,pver),    &! estimate of local volume occupied by clouds
         icwmr1(pcols,pver),  &! in cloud water mixing ration for zhang scheme
         icwmr2(pcols,pver)    ! in cloud water mixing ration for hack  scheme

      real(r8), intent(inout) ::&
         hion(pcols,pver)      ! hydrogen ion concentration (mol/L)

      real(r8), intent(in) ::&
         ekh2o2(pcols,pver)    ! H2O2  Henry's Law coefficient
      real(r8), intent(inout) ::&
         ekso2(pcols,pver)     ! SO2  effective Henry's Law coefficient

      real(r8), intent(out) ::&
         so2new(pcols,pver),  &! so2 array
         so4new(pcols,pver),  &! so4 from aqueous chemistry 
         h2o2new(pcols,pver), &! h2o2 from aqueous chemistry 
         so4srca(pcols,pver), &! so4 production rate from achem
         h2o2snka(pcols,pver)  ! aqueous sink of H2O2

      ! Local variables:

       real(r8)         &
         theff(pcols),  & ! ekso2 at cloudy grid points
         sheff            ! temp array for theff

      integer :: i, k, n
      integer :: i2         ! indice for cloudy grid points
      integer :: ilw        ! number of cloudy grid points
      integer, parameter :: ncyc = 20 ! number of times to loop through aqchem
                                      ! perform aqchem on smaller dt

      real(r8), dimension(pcols,pver) ::&
         cwat,       &! cloud water
         rwat,       &! rain
         twat         ! total liquid water = cwat + rwat
      real(r8) ::&
         h2o2rate,   &! rate of S(IV) + H2O2 reaction 
         o3rate,     &! first order rate coefficient for S(IV) + O3
         molarso4,   &! total so4 in units mol SO4/L H2O
         ahso3        ! (molar)**2 concentration of HSO3-

      real(r8), dimension(pcols) ::&
         tso2,       &! so2 at cloudy grid points
         tso4,       &! so4 at cloudy grid points
         th2o2,      &! h2o2 at cloudy grid points
         to3,        &! o3 at cloudy grid points
         ttwat,      &! twat at cloudy grid points
         temp,       &! tair at cloudy grid points
         tprs,       &! pres at cloudy grid points
         tso2n,      &! so2new at cloudy grid points
         tso4n,      &! so4new at cloudy grid points
         th2o2n,     &! h2o2new at cloudy grid points
         tso4src,    &! so4srca at cloudy grid points
         th2o2snk,   &! h2o2snk at cloudy grid points
         thion,      &! hion at cloudy grid points
         tkh2o2       ! ekh2o2 at cloudy grid points

      real(r8) ::&
         molso2,     &! total so2 in units mol SO2/mol air
         molh2o2,    &! total h2o2 in units mol H2O2/mol air
         molarh2o2,  &! total h2o2 in units mol H2O2/L H2O
         molo3,      &! total o3 in units mol O3/mol air
         eko3,       &! O3   Henry's Law coefficient
         so4src,     &! source of SO4 from SO2 + H2O2
         so2snk,     &! sink of SO2 from SO2 + H2O2
         h2o2snk,    &! sink of H2O2 from SO2 + H2O2
         siv,        &! concentration of S(IV) (mol/L)
         a1, a2,     &! fraction of S(IV) that is HSO3, SO3
         totsb,      &! sum of sulfur at beg of aqchem for one point
         totse,      &! sum of sulfur at end of aqchem for one point
         weight,     &
         tc           ! temperature in deg C

      real(r8) ::&
         dtchem,     &! chemistry timestep
         patm,       &! pressure (atm)
         fa1,        &!  = Khs*K1s/[H+]
         fa2,        &!  = Khs*K1s*K2s/[H+]**2
         khs,        &! Henry's Law constant for SO2
         khsk1s,     &! Khs * K1s
         kh12s,      &! Khs * K1s * K2s
         kh2o2,      &! rate constant for H2O2 + HSO3-
         ko3hs,      &! rate constant for O3 + HSO3-
         ko3so3,     &! rate constant for O3 + SO3=
         qso2,       &! temp array for so2
         qso4,       &! temp array for so4
         qh2o2,      &! temp array for h2o2
         ahion,      &! temp array for ahion
         molrate,    &! rate of rxn in mol/mol-s
         tmprate,    &! temp rate of rxn
         sumrate,    &! h2o2rate + o3rate
         omsm

      ! Function:
!-drb       real(r8) e298, dhr, tt, ek
      real(r8) e298, dhr, tt
      real (r8) ek
      ek(e298, dhr, tt) = e298*exp(dhr*(1._r8/tt - 1._r8/298._r8))
      !----------------------------------------------------------------------

!      call t_startf ('aqchem')

      dtchem = ztodt/ncyc
      omsm = 1._r8 - 1.e-10_r8

      do k=1,pver
         do i=1,ncol

            ! Determine cloud water and rain water mixing ratios
            tc     = tair(i,k) - tmelt
            weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice

            cwat(i,k) = (qc(i,k)/max(cldf(i,k), 0.01_r8) + &
            ! add suspended water from convection and do aqchem in only the liquid phase
                        icwmr1(i,k) + icwmr2(i,k)) * (1._r8-weight)

            if(tair(i,k) .gt. tmelt) then
               weight = 0._r8                    ! fraction of condensate that is ice
            else
               weight = 1._r8
            endif

            rwat(i,k) = qr(i,k)/max(cldv(i,k), 0.01_r8) &
                        *(1._r8-weight) ! aqchem in only the liquid phase

            ! Sum cwat and rwat 
            twat(i,k) = cwat(i,k) + rwat(i,k)

            !        if(twat(i,k) .gt. 0.5e-3) then
            !         tc     = tair(i,k) - tmelt
            !         weight = max(0._r8,min(-tc*0.05,1.0_r8)) 
            !         write(iulog,'(a,3i4,5e12.5)') 'TWAT',i,lat,k, twat(i,k), cwat(i,k),
            !     _    rwat(i,k), icwmr1(i,k)*(1.-weight), icwmr2(i,k)*(1.-weight)
            !         write(iulog,*) ' rwat stuff, qr, cldv, cldf ', qr(i,k), cldv(i,k),
            !     $        cldf(i,k)
            !        endif

         end do
      end do

      ! Initialize arrays
      do k=1,pver
         do i=1,ncol
            so2new(i,k) = so2(i,k)       !keep concentration
            so4new(i,k) = so4(i,k)       
            h2o2new(i,k) = h2o2(i,k)       
            so4srca(i,k) = 0._r8
            h2o2snka(i,k) = 0._r8
         end do
      end do

      ! Initialize cloudy grid point arrays
      !      write(iulog,*) ' Beginning of aqchem'
      do k=1,pver
         ilw = ncpts(k)

         do i2=1,ilw
            tso2(i2)   = so2(ind(i2,k),k)
            tso4(i2)   = so4(ind(i2,k),k)
            th2o2(i2)  = h2o2(ind(i2,k),k)
            to3(i2)    = o3(ind(i2,k),k)
            ttwat(i2)  = twat(ind(i2,k),k)
            temp(i2)   = tair(ind(i2,k),k)
            tprs(i2)   = pres(ind(i2,k),k)
            thion(i2)  = hion(ind(i2,k),k)
            theff(i2)  = ekso2(ind(i2,k),k)
            tkh2o2(i2) = ekh2o2(ind(i2,k),k)
            tso4src(i2) = 0._r8
            th2o2snk(i2) = 0._r8
            !        write(iulog,'(3i5,4e12.5,f8.3,a)') ind(i2,k), lat, k, ttwat(i2), &
            !        thion(i2), tso4(i2), tso2(i2), -log10(thion(i2)), '  begin'
         end do
                    
         ! Aqueous chemistry calculations begin here
         do i2=1,ilw
            patm = tprs(i2)/101325._r8
            khs    = ek(1.23_r8,3120._r8,temp(i2))
            khsk1s = khs * ek(1.3e-2_r8,2015._r8,temp(i2))
            kh12s = khsk1s * ek(6.31e-8_r8,1505._r8,temp(i2)) 
            kh2o2  = 6.2534e14_r8*exp(-4751._r8/temp(i2)) 
            molo3 = to3(i2) 
            eko3  = ek(1.15e-2_r8,2560._r8,temp(i2))
            ko3hs = 4.28e13_r8*exp(-5533._r8/temp(i2)) 
            ko3so3= 7.436e16_r8*exp(-5280._r8/temp(i2)) 

            !        if(lat .eq. 8 .and. k .eq. 14 .and. i2 .gt. 8) then
            !         write(iulog,'(i3,7e11.4)') i2, tso2(i2), th2o2(i2), thion(i2), &
            !          theff(i2), temp(i2), patm, molo3
            !         write(iulog,'(3x,e11.4,f8.1)') ttwat(i2), dtchem
            !        endif
!cdir expand=ncyc
            do n=1,ncyc                     !Do aqchem on a smaller dt   5/5/97

               if(n .eq. 1) then
                  qso2 = tso2(i2)
                  qso4 = tso4(i2)
                  qh2o2 = th2o2(i2)
                  ahion = thion(i2)
                  sheff = theff(i2)
               endif

               fa1 = khsk1s/max(ahion,1.e-30_r8)
               fa2 = kh12s/max((ahion*ahion),1.e-30_r8)

               ! Calculate rate of HSO3 + H2O2 reaction
               molso2  = qso2  * mwa2so2
               molh2o2 = qh2o2 * mwa2h2o2
               ahso3 = fa1 * patm*molso2
               molarh2o2= tkh2o2(i2)*patm*molh2o2

               h2o2rate = kh2o2 * ahion * ahso3 * molarh2o2/(1._r8 + 13._r8*ahion)  

               ! Convert to mol/(mol-s) 
               molrate = h2o2rate*ttwat(i2)*28.97e-3_r8
               ! check for conservation in kg/(kg-s)
               tmprate = mwa2h2o2*min(molrate/mwa2h2o2, qh2o2/dtchem)
               tmprate = mwa2so2*min(tmprate/mwa2so2, qso2/dtchem)

               h2o2rate = tmprate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3_r8)

               ! Calculate rate of o3 oxidation
               siv = sheff * molso2 * patm         ![S(IV)]
               a1 = fa1 / sheff
               a2 = fa2 / sheff

               o3rate  = (ko3hs*a1 + ko3so3*a2) * eko3*molo3 * patm*siv

               ! Convert to mol/(mol-s) 
               molrate = o3rate*ttwat(i2)*28.97e-3_r8
               ! check for conservation in kg/(kg-s)
               tmprate = mwa2so2*min(molrate/mwa2so2, qso2/dtchem)

               o3rate = tmprate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3_r8)

               ! check total rate for conservation
               sumrate = h2o2rate + o3rate
               !         if(lat .eq. 8 .and. k .eq. 14 .and. i2 .gt. 8) then
               !          write(iulog,'(i3,5e11.4)') n, sumrate, h2o2rate, o3rate, &
               !            qso2/dtchem, qh2o2/dtchem
               !         endif
               molrate = sumrate*ttwat(i2)*28.97e-3_r8
               molrate = mwa2so2*min(molrate/mwa2so2, qso2/dtchem)
               tmprate = molrate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3_r8)
               h2o2rate = h2o2rate/max(sumrate, 1.e-30_r8) * tmprate
               o3rate = o3rate/max(sumrate, 1.e-30_r8) * tmprate

               ! Update so2, so4, h2o2
               so2snk = (h2o2rate + o3rate)*ttwat(i2)*1.e-3_r8 * 64._r8
               h2o2snk = h2o2rate*ttwat(i2)*1.e-3_r8 * 34._r8
               !         if(qh2o2 .lt. omsm*h2o2snk*dtchem) then
               !          write(iulog,'(a,2e12.5,3i4)') 'AQCHEM qh2o2 < h2o2snk*dt ', qh2o2, &
               !           h2o2snk*dtchem, lat, ind(i2,k), k
               !         endif
               !         if(qso2 .lt. omsm*so2snk*dtchem) then
               !          write(iulog,'(a,2e12.5,3i4)') 'AQCHEM qso2 < so2snk*dt ', qso2, &
               !           so2snk*dtchem, lat, ind(i2,k), k
               !         endif
               so4src = 1.5_r8*so2snk

               qso2 = max(qso2 - so2snk*dtchem, 0._r8)
               qso4 = qso4 + so4src*dtchem
               qh2o2 = max(qh2o2 - h2o2snk*dtchem, 0._r8)

               tso4src(i2) = tso4src(i2) + so4src
               th2o2snk(i2) = th2o2snk(i2) + h2o2snk

               ! calculate new ahion, sheff
               molso2 = qso2 * mwa2so2
               molarso4 = qso4/(max(ttwat(i2), 1.e-30_r8)*96.e-3_r8) 
               ahso3 = khsk1s * molso2 * patm
!++pjr
! rewrote this expression so we dont take the sqrt(max(a,1.e-60))
!                               but rather max(sqrt(a),1.e-30)
!               ahion = (molarso4 + sqrt(max(molarso4*molarso4 + 4.*ahso3, &
!                       1.e-60_r8)))*0.5
!--pjr
               ahion = (molarso4 + sqrt(max(molarso4*molarso4 + 4._r8*ahso3, &
                       0._r8)))*0.5_r8
               ! Limit to keep pH > 0.0 (ion mixing ratio < 1)
               ahion = min(max(ahion,1.e-30_r8),0.999_r8)

               sheff = khs + khsk1s/ahion + kh12s/max((ahion*ahion),1.e-30_r8)
            end do    !ncyc

            tso2n(i2) = qso2
            tso4n(i2) = qso4
            th2o2n(i2)= qh2o2
            thion(i2) = ahion
            theff(i2) = sheff
            tso4src(i2) = tso4src(i2)/ncyc
            th2o2snk(i2) = th2o2snk(i2)/ncyc

         end do


         do i2=1,ilw
            ! so2new is modified for just the cloudy portion of the grid cell.  Now
            !  adjust it for the entire grid cell.
            so2new(ind(i2,k),k) = cldv(ind(i2,k),k)*tso2n(i2) + &
                                  (1._r8-cldv(ind(i2,k),k))*tso2(i2)
            so4new(ind(i2,k),k) = cldv(ind(i2,k),k)*tso4n(i2) + &
                                  (1._r8-cldv(ind(i2,k),k))*tso4(i2)
            so4srca(ind(i2,k),k) = cldv(ind(i2,k),k)*tso4src(i2) 
            h2o2snka(ind(i2,k),k)= cldv(ind(i2,k),k)*th2o2snk(i2) 
            h2o2new(ind(i2,k),k) = cldv(ind(i2,k),k)*th2o2n(i2) + &
                                   (1._r8-cldv(ind(i2,k),k))*th2o2(i2)
            hion(ind(i2,k),k) = thion(i2)  
            ekso2(ind(i2,k),k) = theff(i2)  
            !         write(iulog,'(3i5,4e12.5,f8.3,a)') ind(i2,k), lat, k, ttwat(i2), &
            !         thion(i2), tso4(i2), tso2(i2), -log10(thion(i2)), '  end'
         end do


         do i=1,ncol
            totsb = 0.5_r8*so2(i,k) + 1._r8/3._r8 * so4(i,k)
            totse = 0.5_r8*so2new(i,k) + 1._r8/3._r8 * so4new(i,k)
            !        if(abs(totsb-totse) .gt. .001*totsb) then
            !         write(iulog,*) 'Total sulfur not conserved in aqueous chemistry'
            !         write(iulog,*) i, k, totsb, totse
            !         stop
            !        endif

         end do

      end do

!      call t_stopf ('aqchem')
      return
   end subroutine aqchem

!##############################################################################

   subroutine cldychmini( tair, pres, qc, qr, cldf, &
                          cldv, icwmr1, icwmr2, so2, so4, &
                          ph, hion, ekso2, ekh2o2, ind, &
                          ncpts, ncol ) 

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Calculate the pH, solubility constants of SO2 and H2O2, and the
      !  location of the cloudy grid points
      ! 
      ! Author: M. Barth
      !----------------------------------------------------------------------- 

      implicit none

      real(r8), intent(in) ::&
         tair(pcols,pver),    &! air temperature (K)
         pres(pcols,pver),    &! midpoint pressures
         qc(pcols,pver),      &! cloud water mixing ratio
         qr(pcols,pver),      &! rain mixing ratio
         cldf(pcols,pver),   &! fraction cloud cover
         cldv(pcols,pver),    &! cloudy volume which is occupied by rain or cloud water
         so2(pcols,pver),     &! so2 array
         so4(pcols,pver),     &! so4 array
         icwmr1(pcols,pver),  &!  in cloud water mixing ration for zhang scheme
         icwmr2(pcols,pver)    !  in cloud water mixing ration for hack  scheme

      real(r8), intent(out) ::&
         ph(pcols,pver),      &! pH of drops
         hion(pcols,pver)      ! hydrogen ion concentration (mol/L)

      real(r8), intent(out) ::         &
         ekh2o2(pcols,pver),  &! H2O2 Henry's Law coefficient
         ekso2(pcols,pver)     ! SO2  effective Henry's Law coefficient

      integer, intent(in) ::  &
         ncol                  ! number of atmospheric columns used in chunk
      integer, intent(out) ::&
         ind(pcols,pver),     &! indices for cloudy grid points
         ncpts(pver)           ! number of cloudy grid points

      ! Local variables:

      real(r8) :: &
         theff(pcols)          ! ekso2 at cloudy grid points

      integer :: &
         i, k,                &
         i2,                  &! index for cloudy grid points
         ilw                   ! number of cloudy grid points

      real(r8), dimension(pcols,pver) ::&
         cwat,       &! cloud water
         rwat,       &! rain
         twat         ! total liquid water = cwat + rwat

      real(r8), dimension(pcols) ::&
         tso2,       &! so2 at cloudy grid points
         tso4,       &! so4 at cloudy grid points
         ttwat,      &! twat at cloudy grid points
         temp,       &! tair at cloudy grid points
         tprs,       &! pres at cloudy grid points
         tph,        &! ph at cloudy grid points
         thion        ! hion at cloudy grid points
      real(r8)  ::   &
         molarso4,   &! total so4 in units mol SO4/L H2O
         ahso3,      &! (molar)**2 concentration of HSO3-
         molso2,     &! so2 in mol/mol units
         weight,     &
         tc           ! temperature in deg C

      ! Function:
!-drb      real(r8) e298, dhr, tt, ek
      real(r8) e298, dhr, tt
      real(r8) ek
      ek(e298, dhr, tt) = e298*exp(dhr*(1._r8/tt - 1._r8/298._r8))
      !----------------------------------------------------------------------

      do k=1,pver
         do i=1,ncol

            ekh2o2(i,k) = ek(7.4e4_r8, 6621._r8, tair(i,k))
            ekso2(i,k)  = ek(1.23_r8,  3120._r8, tair(i,k)) 
            ! Determine cloud water and rain water mixing ratios
            tc     = tair(i,k) - tmelt
            weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice

            cwat(i,k) = (qc(i,k)/max(cldf(i,k), 0.00001_r8) + &
               ! add suspended water from convection and do aqchem in only the liquid phase
                        icwmr1(i,k) + icwmr2(i,k)) * (1._r8-weight)

            if(tair(i,k) .gt. tmelt) then
               weight = 0._r8                    ! fraction of condensate that is ice
            else
               weight = 1._r8
            endif

            rwat(i,k) = qr(i,k)/max(cldv(i,k), 0.00001_r8) &
                        *(1._r8-weight) ! aqchem in only the liquid phase

            ! Sum cwat and rwat 
            twat(i,k) = cwat(i,k) + rwat(i,k)

         end do
      end do

      ! Initialize arrays
      do k=1,pver
         do i=1,ncol
            ph(i,k) = 99._r8
            hion(i,k) = 0._r8
            ind(i,k) = 0
         end do
      end do

      ! Find which grid points have liquid water
      do k=1,pver
         ncpts(k) = 0
         ilw = 0
         do i=1,ncol
            if(cldv(i,k) .ge. 0.02_r8 .and. twat(i,k) .ge. 1.e-12_r8) then
               ilw = ilw + 1
               ind(ilw,k) = i
            endif
         end do

         ncpts(k) = ilw       !assign number of cloudy grid points to array

         do i2=1,ilw
            tso2(i2) = so2(ind(i2,k),k)
            tso4(i2) = so4(ind(i2,k),k)
            ttwat(i2) = twat(ind(i2,k),k)
            temp(i2) = tair(ind(i2,k),k)
            tprs(i2) = pres(ind(i2,k),k)

            ! Set pH as a function of HSO3- and SO4=
            !  Say NH4+ = 1.5 SO4= and NO3- = 0.5 SO4=
            !  Thus, H+ = HSO3- + SO4=  as done in ECHAM
            molso2 = tso2(i2)  * mwa2so2

            molarso4 = tso4(i2)/(max(ttwat(i2), 1.e-30_r8)*96.e-3_r8) 
            ahso3 = ek(1.23_r8,3120._r8,temp(i2)) * &
                    ek(1.3e-2_r8,2015._r8,temp(i2)) * &
                    molso2 * tprs(i2)/101325._r8
!++pjr
! rewrote this expression so we dont take the sqrt(max(a,1.e-50))
!                               but rather max(sqrt(a),1.e-25)
!            thion(i2) = (molarso4 + sqrt(max(molarso4*molarso4 + 4.*ahso3, &
!                         1.e-50_r8)))*0.5
!--pjr
            thion(i2) = (molarso4 + &
                         sqrt(max(molarso4*molarso4 + 4._r8*ahso3,0._r8)))*0.5_r8


!++ag Limit to keep pH > 0.0 (ion mixing ratio < 1)
!           thion(i2) = max(thion(i2),1.e-25_r8)
            thion(i2) = min(max(thion(i2),1.e-25_r8),0.999_r8)
!--ag
            tph(i2) = -log10(thion(i2))

            theff(i2) = ek(1.23_r8, 3120._r8, temp(i2)) * (1._r8 + &
                        ek(1.3e-2_r8, 2015._r8, temp(i2))/thion(i2) * (1._r8 + &
                        ek(6.31e-8_r8, 1505._r8, temp(i2))/thion(i2)))

         end do


         do i2=1,ilw
            hion(ind(i2,k),k) = thion(i2)  
            ph(ind(i2,k),k) = tph(i2)  
            ekso2(ind(i2,k),k) = theff(i2)

         end do
      end do

   end subroutine cldychmini


!##############################################################################

   subroutine coszen( calday, dodiavg, clat, clon, ncol, coszrs )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Compute solar distance factor and cosine solar zenith angle usi
      ! day value where a round day (such as 213.0) refers to 0z at
      ! Greenwich longitude.
      !
      ! Use formulas from Paltridge, G.W. and C.M.R. Platt 1976: Radiative
      ! Processes in Meterology and Climatology, Elsevier Scientific
      ! Publishing Company, New York  p. 57, p. 62,63.
      ! 
      ! Author: CCM2
      !----------------------------------------------------------------------- 

      implicit none

      real(r8), intent(in) ::&
         calday,              &! Calendar day, including fraction
         clat(pcols),         &! latitudes (radians)
         clon(pcols)

      integer, intent(in) ::&
!         lat(pcols),          &! latitude indices
!         lon(pcols),          &! longitude indices
         ncol                  ! number of atmospheric columns used in chunk

      logical, intent(in) ::&
         dodiavg               ! true => do diurnal averaging

      real(r8), intent(out) ::&
         coszrs(pcols)        ! Cosine solar zenith angle

      ! Local variables
      integer :: i     ! Longitude loop index
      real(r8) ::&
         phi,     &! Greenwich calendar day + local time + long offset
         theta,   &! Earth orbit seasonal angle in radians
         delta,   &! Solar declination angle  in radians
         sinc,    &! Sine   of latitude
         cosc,    &! Cosine of latitude
         sind,    &! Sine   of declination
         cosd      ! Cosine of declination
      real(r8) ::&
         frac,    &! Daylight fraction
         arg,     &!
         tsun,    &! temporary term in diurnal averaging
         coszrsu   ! uniform cosine zenith solar angle 
      !-----------------------------------------------------------------------

      theta = 2._r8*pie*calday/dayspy

      ! Solar declination in radians:
      delta = .006918_r8 - .399912_r8*cos(theta) + .070257_r8*sin(theta) - &
              .006758_r8*cos(2._r8*theta) + .000907_r8*sin(2._r8*theta) - &
              .002697_r8*cos(3._r8*theta) + .001480_r8*sin(3._r8*theta)

      do i=1,ncol
         ! Compute local cosine solar zenith angle,
         sinc = sin(clat(i))
         sind = sin(delta)
         cosc = cos(clat(i))
         cosd = cos(delta)
   
         ! If using diurnal averaging, then compute the average local cosine solar 
         ! zenith angle using formulas from paltridge and platt 1976  p. 57, p. 62,63.
         if (dodiavg) then
            arg = -(sinc/cosc)*(sind/cosd)
            if (arg .lt. -1._r8) then
               frac = 1.0_r8
            else if (arg .gt. 1._r8) then
               frac = 0.0_r8
            else
               frac = (1._r8/pie)*acos(arg)
            endif
            tsun = pie*frac
            if (tsun .gt. 0._r8) then
               coszrsu =  sinc*sind + (cosc*cosd*sin(tsun))/tsun
            else
               coszrsu = 0.0_r8
            endif
            coszrs(i) = coszrsu
         else                       ! No diurnal averaging
   
            ! Calday is the calender day for Greenwich, including fraction
            ! of day; the fraction of the day represents a local time at
            ! Greenwich; to adjust this to produce a true instantaneous time
            ! For other longitudes, we must correct for the local time change:
            ! local time based on the longitude and day of year
            ! then compute the local cosine solar zenith angle
	     coszrs(i) = sinc*sind - cosc*cosd*cos(2._r8*pie*calday + clon(i))
         end if
      end do

   end subroutine coszen

!##############################################################################

   subroutine dicor( calday, clat, ncol, corr )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Calculate a correction factor for the ho2 reaction rate to account for
      ! the diurnal variation of ho2
      !
      ! Author: Mary Barth
      !----------------------------------------------------------------------- 

      implicit none

      real(r8), intent(in) :: calday                ! Calendar day, including fraction
      real(r8), intent(in) :: clat(pcols)           ! latitudes (radians)

      integer, intent(in) ::  ncol                  ! number of atmospheric columns used in chunk

      real(r8), intent(out) :: corr(pcols)          ! correction factor

      !---------------------------Local variables-----------------------------

      integer ic        ! index in chunk
      integer i         ! Longitude loop index
      real(r8) phi      ! Greenwich calendar day + local time + long offset
      real(r8) theta    ! Earth orbit seasonal angle in radians
      real(r8) delta    ! Solar declination angle  in radians
      real(r8) sinc     ! Sine   of latitude
      real(r8) cosc     ! Cosine of latitude
      real(r8) sind     ! Sine   of declination
      real(r8) cosd     ! Cosine of declination
      real(r8) coszrs       ! Cosine solar zenith angle
      real(r8) xxx      ! Work space
      real(r8) yyy      ! Work space
      integer nc        ! counter

      real(r8) ::     dicor_corr          ! h2o reaction rate correction
      integer ntimes    ! number of points to evaluate the diurnal average

      !-----------------------------------------------------------------------
      !

      do ic=1,ncol

         ! Compute solar distance factor and cosine solar zenith angle using
         ! day value where a round day (such as 213.0) refers to 0z at
         ! Greenwich longitude.
         !
         ! Use formulas from Paltridge, G.W. and C.M.R. Platt 1976: Radiative
         ! Processes in Meterology and Climatology, Elsevier Scientific
         ! Publishing Company, New York  p. 57, p. 62,63.

         theta = 2._r8*pie*calday/dayspy
         !
         ! Solar declination in radians:
         !
         delta = .006918_r8 - .399912_r8*cos(theta) + .070257_r8*sin(theta) - &
              .006758_r8*cos(2._r8*theta) + .000907_r8*sin(2._r8*theta) - &
              .002697_r8*cos(3._r8*theta) + .001480_r8*sin(3._r8*theta)

         ! now calculate the solar zenith angle 
         ! and some useful quantities for the whole day
         !
         ! Compute local cosine solar zenith angle,
         !
         sinc = sin(clat(ic))
         sind = sin(delta)
         cosc = cos(clat(ic))
         cosd = cos(delta)
         !
         ! Calday is the calender day for Greenwich, including fraction
         ! of day; 
         ! since all we care about in this calculation is the diurnal variation
         ! of some quantities we just increment calday over one day.
         !
         xxx = 0._r8
         yyy = 0._r8
         nc  = 0
         ntimes = 48
         do i=1,ntimes
            phi       = calday + (real(i,r8)-1)/real(ntimes,r8)
            coszrs = sinc*sind - cosc*cosd*cos(2._r8*pie*phi)
            if (coszrs.gt.0) then
               nc  = nc + 1
               xxx = xxx + coszrs**2
               yyy = yyy + coszrs
            endif
         end do

         if (yyy.gt.0._r8) then
            dicor_corr = ntimes*xxx/(yyy**2)
         else
            dicor_corr = 1._r8
         endif


         !
         ! Return correction
         !
         corr(ic) = dicor_corr
      end do

    end subroutine dicor

!##############################################################################

   subroutine gaschem( so2, so4, dms, h2o2, tair, &
                       pres, so2new, so4new, dmsnew, h2o2new, &
                       h2o23d, oh, no3, ho2, h2o, &
                       jh2o2, so4srcg, dmssnk, so2srcg, so2srcg2, &
                       h2o2src, h2o2snkg, ztodt, ncol )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Calculate the sources and sinks from the gas chemistry:
      !          DMS + OH     --> a*SO2 + (1-a)*MSA
      !          DMS + NO3    --> SO2 
      !          SO2 + OH + M --> SO4 + M
      !          HO2 + HO2    --> H2O2    
      !          H2O2+ OH     --> HO2 + H2O   
      !          H2O2+ hv     --> 2OH
      !
      !  where a is the yield of SO2 from DMS oxidation, MSA is methane sulfonic
      !   acid which is not carried in this model, and M is an air molecule and
      !   is accounted for in the rate coefficient expression.
      !
      ! Author: M. Barth
      !----------------------------------------------------------------------- 

      implicit none

      real(r8), intent(in) :: &
         so2(pcols,pver),     &! so2 array
         so4(pcols,pver),     &! so4 from gas chemistry 
         dms(pcols,pver),     &! dms
         h2o2(pcols,pver),    &! h2o2
         oh(pcols,pver),      &! oh
         no3(pcols,pver),     &! no3
         ho2(pcols,pver),     &! ho2
         h2o(pcols,pver),     &! h2o
         jh2o2(pcols,pver),   &! photolysis rate of H2O2
         tair(pcols,pver),    &! air temperature (K)
         pres(pcols,pver),    &! midpoint pressures
         ztodt,               &! time step
         h2o23d(pcols,pver)    ! h2o2 -- 3d field from images

      integer, intent(in) ::  &
!         lat(pcols),          &! latitude indices
         ncol                  ! number of atmospheric columns used in chunk

      real(r8), intent(out) ::  &
         so4new(pcols,pver),    &! so4 from gas chemistry 
         dmsnew(pcols,pver),    &! dms
         h2o2new(pcols,pver),   &! h2o2
         so4srcg(pcols,pver),   &
         so2srcg2(pcols,pver)

      real(r8), intent(inout) ::  &
         so2new(pcols,pver),      &! so2 array
         dmssnk(pcols,pver),      &
         so2srcg(pcols,pver),     &
         h2o2src(pcols,pver),     &
         h2o2snkg(pcols,pver)

      ! Local variables:

      integer ::&
         i, k, ik
      logical :: found
      real(r8), parameter ::&         !  from NASA(1992)
         rl300=3.e-31_r8,    &! rate constant at low pressure
         an=3.3_r8,          &! constant:  (300/Tair)**an
         rh300=1.5e-12_r8,   &! rate constant at high pressure
         am=0._r8             ! constant:  (300/Tair)**am

      real(r8) ::&
         air,             &! concentration of air (molec/cm3)
         rlt,             &! used in calculation of SO2 rate coefficient
         rht,             &! used in calculation of SO2 rate coefficient
         ax, bx,          &! used in calculation of SO2 rate coefficient
         psi,             &! used in calculation of SO2 rate coefficient
         so2rate,         &! first order rate coefficient for SO2 + OH
         so2snk,          &
         dmsrate1,        &! first order rate coefficient for DMS + OH
         dmsno3,          &! first order rate coefficient for DMS + NO3
         dmstot            ! first order rate coefficient for both rxns
      real(r8) rk2         ! rate coeff. in HO2 + HO2 rxn
      real(r8) rk3         ! rate coeff. in HO2 + HO2 rxn
      real(r8) fh2o        ! rate coeff. in HO2 + HO2 rxn
      real(r8) ho2rate     ! total rate coeff for HO2 + HO2 rxn
      real(r8) h2o2des     ! rate for H2O2*OH
      real(r8) h2o2phot    ! rate jh2o2*H2O2
      !      real(r8) tauh2o2         !relaxation time for H2O2 generation
      !      real(r8) xh2o2           !zonally avgd H2O2 in kg/kg
      !----------------------------------------------------------------------c
      
      !  Assume a 1.5 day relaxation time for H2O2 generation
      !      tauh2o2 = 129600.    !seconds

!      call t_startf ('gaschem')

      found = .false.
      do k=1,pver
         do i=1,ncol
            ! SO2 + OH rate coefficient:
            if(tair(i,k) .le. 0._r8) then
               found = .true.
               exit
            endif
         end do
      end do
      if (found) then
        do k=1,pver
           do i=1,ncol
              ! SO2 + OH rate coefficient:
              if(tair(i,k) .le. 0._r8) then
                 write(iulog,*) 'TAIR bad ', i, k, tair(i,k)
                 call endrun ('GASCHEM')
              endif
           end do
        end do
      endif
      
      do k=1,pver
         do i=1,ncol

            so2rate = 0._r8

            ! determine concentration of air (molec/cm3) = rho(kg/m3)*Na/MWair
            air = pres(i,k)/(287._r8*tair(i,k)) * 1.e-3_r8/28.966_r8 * 6.022e23_r8
            !         air = 0.8*air                   !N2 concentration
 
            ! SO2 + OH rate coefficient:
            rlt = rl300 * (300._r8/tair(i,k))**an
            rht = rh300 * (300._r8/tair(i,k))**am
            ax = rlt*air/rht 
            psi = 1._r8/(1._r8 + (log10(ax))**2)
            bx = rlt*air/(1._r8 + ax)
            so2rate = bx * 0.6_r8**psi             !rate coef. in cm3/molec-s
            ! rewrite so2rate as first order rate coefficient
            so2rate = so2rate * oh(i,k)

            so2snk = so2(i,k)*(1._r8-exp(-so2rate*ztodt))
            ! Update sulfur species
            so2new(i,k) = so2(i,k)  - so2snk

            so4new(i,k) = so4(i,k) + 1.5_r8*so2snk
            so4srcg(i,k) = 1.5_r8*so2snk/ztodt

            !   Note: 1.5 factor accounts for different molecular weights of so4 and so2
            !    mw(so4)=96;  mw(so2)=64;  mw(so4)/mw(so2)=3/2=1.5
         end do
      end do

      do k=1,pver
         do i=1,ncol

            ! DMS oxidation
            !c DMS + OH rate coefficient and SO2 yield:
            dmstot = 0._r8
            dmsrate1 = 0._r8
            dmsno3 = 0._r8
            ! Get dmsrate from tabulated data
            ik = int(tair(i,k)-195._r8)
            ik = min(max(ik,1), 121)

            if(dms(i,k) .ge. 1.e-30_r8) then
               dmsrate1 = dmsrate(ik)*oh(i,k)

               ! DMS + NO3
               dmsno3 = 1.e-12_r8*exp(-500._r8*(0.0033557_r8 - 1._r8/tair(i,k))) * &
                        no3(i,k)

               ! Total first order rate coeff for DMS
               dmstot = dmsrate1 + dmsno3
            endif

            ! Update sulfur species
            dmssnk(i,k) = dms(i,k)*(1._r8-exp(-dmstot*ztodt)) / ztodt
            dmsnew(i,k) = dms(i,k) - dmssnk(i,k)*ztodt
            so2srcg(i,k) = 1.032258_r8 * (yield(ik) * dms(i,k)*(1._r8-exp(-dmsrate1*ztodt)) + &
                           dms(i,k)*(1._r8-exp(-dmsno3*ztodt)) ) / ztodt
            so2new(i,k) = so2new(i,k) + so2srcg(i,k)*ztodt
            so2srcg2(i,k) = 1.032258_r8 / ztodt * &
                            yield(ik) * dms(i,k)*(1._r8-exp(-dmsrate1*ztodt)) 
            !   Note: 1.032258 = 64/62 accounts for different MW of so2 and dms
         end do
      end do

      do k=1,pver
         do i=1,ncol
            !-----------------------------------------------------------------------
            !c Generate H2O2 via gas phase processes -- parameterizing chemistry
            !        xh2o2 = h2o23d(i,k) * (287.*tair(i,k)/pres(i,k))*34./6.022e20
            !        h2o2new(i,k) = h2o2(i,k) + ztodt/tauh2o2 * (xh2o2 - h2o2(i,k))
            !
            ! Generate H2O2 via the HO2 + HO2 reaction
            !  HO2 concentrations are prescribed
            !  water vapor concentrations taken from q(i,k,1)
            !
            ! determine concentration of air (molec/cm3) = rho(kg/m3)*Na/MWair
            air = pres(i,k)/(287._r8*tair(i,k)) * 1.e-3_r8/28.966_r8 * 6.022e23_r8
            rk2 = 1.7e-12_r8*exp( 600._r8*(1._r8/tair(i,k) - 0.0033557_r8))
            rk3 = 4.9e-32_r8*exp(1000._r8*(1._r8/tair(i,k) - 0.0033557_r8))
            fh2o = 1.4e-21_r8*exp(2200._r8/tair(i,k))
            ho2rate = (rk2 + rk3*air)*(1._r8 + fh2o*h2o(i,k)) * &
                      ho2(i,k)*ho2(i,k)
            !        xh2o2 = h2o23d(i,k) * (287.*tair(i,k)/pres(i,k))*34./6.022e20

            ! Gas-phase destruction of H2O2 occurs via:
            !      H2O2 + hv --> 2OH
            !      H2O2 + OH --> HO2 + H2O
            !  These reactions seem to have the same reaction rate.  Therefore
            !   parameterize this chemistry as doubling the rate of reaction of
            !   H2O2 + OH.
      
            ! h2o2*air*28.966/34. = h2o2 in molec/cm3

            h2o2des = 1.7e-12_r8*exp(-200._r8*(1._r8/tair(i,k) - 0.0033557_r8)) * &
                      oh(i,k) * h2o2(i,k)*air*mwa2h2o2

            h2o2phot = jh2o2(i,k) * h2o2(i,k)*air*mwa2h2o2

            h2o2src(i,k) = ho2rate/(air*mwa2h2o2)
            h2o2snkg(i,k) = (h2o2des + h2o2phot)/(air*mwa2h2o2)
!           h2o2new(i,k) = h2o2(i,k) + (ho2rate - h2o2des - h2o2phot) * &
!                           ztodt/(air*mwa2h2o2)
            h2o2new(i,k) = h2o2(i,k) + (h2o2src(i,k) - h2o2snkg(i,k)) * ztodt

            !c        h2o2new(i,k) = min(h2o2(i,k) + 
            !c     _      ho2rate * ztodt*(287.*tair(i,k)/pres(i,k))*34./6.022e20,
            !c     _      xh2o2)

         end do
      end do

!      call t_stopf ('gaschem')
      
      return
   end subroutine gaschem

!##############################################################################

   subroutine getj( clat, clon, calday, zm, ncol, jout )

      !----------------------------------------------------------------------- 
      ! Purpose: 
      ! Set h2o2 photolysis rates.
      ! 
      ! Author: M. Barth
      !----------------------------------------------------------------------- 

      implicit none

      real(r8), intent(in) ::&
         clat(pcols),           &! latitudes (radians)
         clon(pcols),           &! longitude (radians)
         calday,                &! Current calendar day = julian day + fraction
         zm(pcols,pver)          ! height of midpoint of layer

      integer, intent(in) ::&
!         lat(pcols),            &! latitude indices
!         lon(pcols),            &! longitude indices
         ncol                    ! number of atmospheric columns used in chunk

      real(r8), intent(out) ::&
         jout(pcols,pver)        ! photolysis rate for location

      ! Local variables
      logical dodiavg                 ! when true, diurnally avg zenith angle
      real(r8) cosza(pcols)               ! cosine of zenith angle
      real(r8) dicosza                    ! diurnally avgd cosine of zenith angle
      real(r8) disecza                    ! diurnally avgd cosine of zenith angle

      integer i, is, k
      integer ind(3)
      integer mz
      real(r8) psum
      real(r8) ratza                      !ratio of zenith angles
      real(r8) ratalb                     !ratio of albedos = (0.3-0.2)/(0.5-0.2)
      real(r8) sumrat                     !sum of the ratios
      real(r8) w1                         !1-sumrat
      !-----------------------------------------------------------------------c
      dodiavg = .true.
      call coszen( calday, dodiavg, clat, clon, ncol, cosza )
! TBH:  Loops are not in optimal order for memory access...  
      do i=1,ncol
         dicosza = cosza(i)
   
         if( dicosza .lt. 0.02_r8 ) then
            do k=1,pver
               jout(i,k) = 0._r8
            end do
         else
            disecza = min( 1._r8/max( dicosza, 1.e-7_r8 ), 50._r8 )
      
            ! Interpolate data onto current level and zenith angle for alb=0.3
      
            ! zenith angle:
            do is = 1, 8                            !zasec = za in table
               if( zasec(is) .gt. disecza ) go to 5 !secza = za at gridpt
            end do
            is = is - 1
5           is = max( is-1, 1 )
            ratza = (disecza - zasec(is)) * delang(is)
      
            !     albedo -- set to 0.3 (data is for 0.2 and 0.5)
            ratalb = (0.3_r8 - alb(1)) * delalb
      
            sumrat = ratza + ratalb
            w1 = 1._r8 - sumrat
            
            ! Set indices
            ! note albedo index is either 1 or 2; therefore only one jh2o2 will
            ! have ialb=2
            ind(1) = nz2*(is-1) 
            ind(2) = nz2*is 
            ind(3) = nz2*nza + nz2*(is-1) 
      
            do k = 1, pver
               mz = nint( min( 50._r8, zm(i,k)*0.001_r8 + 1._r8 ) )
               psum = w1*jper2(ind(1)+mz) + ratza*jper2(ind(2)+mz) + &
                      ratalb*jper2(ind(3)+mz)
               jout(i,k) = exp( psum )
            end do
         endif
      end do

   end subroutine getj

!##############################################################################

end module sulchem

